

library(HMM)

len=500
dim1=3
startt=1
endd=3
# offsett=.01


stats=as.character(seq(startt, endd))
symbolss=c('a')
# startP=c(0.01, 0.04, .9, 0.04, .01)
startP=c(.1, .99, .1)
startP=startP/sum(startP)

# mean0=0;
# sigm=.1;


# transP=matrix(rep(0,dim1^2), nrow=dim1)
emissionP=matrix(rep(1/length(symbolss),dim1*length(symbolss)), nrow=dim1)





pdf('1.cn_states_simu.pdf',width=9, height=6)



p11=1-seq(0, .3, .01)
p22=1-seq(0, .1, .003)

obs=rep(symbolss, len)

for (i in 1:length(p11)) {
	for (j in 1:length(p22)) {
	
	trans1=c(p11[i], .1, .01)
	trans1=trans1/sum(trans1)
	trans2=c(.1, p22[j], .1)
	trans2=trans2/sum(trans2)
	trans3=c(.01, .1, p11[i])
	trans3=trans3/sum(trans3)
	transP=matrix(c(trans1, trans2, trans3), nrow=dim1, byrow=T)
	
	hmm1=initHMM(stats, symbolss, startProbs=startP, transProbs=transP, emissionProbs=emissionP)
	# simhmm1=simHMM(hmm1, len)
	path1=viterbi(hmm1, obs)
	path1=as.numeric(path1)

	write.table(path1, file = paste(i*j, 'cn_simu', 'txt', sep='.'), quote =F, row.names=F, col.names=F)
	plot(path1, ylab='CN',xlab='INDEX', main='SIMU CN', col='blue', type='s')

	}
}




# for (i in seq(num1[1]*num1[2])) {

# for (i in seq(dim1)) {
# norm1=rnorm(dim1, mean0, sigm); pnorm(norm1, mean0, sigm)
# norm1=pnorm(norm1+offsett)-pnorm(norm1)
# norm1=norm1/sum(norm1)
# transP[i,]=norm1
# }

# for (j in seq(dim1)) {
# norm1=c(.9)


# norm1=rnorm(dim1, mean0, sigm); pnorm(norm1, mean0, sigm)
# norm1=pnorm(norm1+offsett)-pnorm(norm1)
# norm1=norm1/sum(norm1)
# transP[i,]=norm1
# }




dev.off()
